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I. COMPUTATION AND ANALYSIS OF FREE SURFACE FLOW OF A THIN 
LIQUID FILM UNDER ZERO AND NORMAL GRAVITY 


1 . 1 Summary 

The results of a numerical computation and theoretical analysis are 
presented for the flow of a thin liquid film in the presence or absence of 
a gravitational body force. The numerical computations employed a 
curvilinear body- fitted coordinate system. The governing transport 
equations were discretized using the finite control- volume formulation and 
solved implicitly. The computer code PHOENICS was used for this purpose. 

Five different flow systems were studied: (a) falling film down a 
vertical wall, (b) plane film flow along a horizontal plate, (c) plane film 
flow at zero gravity, (d) radial film flow along a horizontal plate and (e) 
radial film flow at zero gravity. In all these systems, the distribution 
of the film height, the flow field and resistance exerted by the solid wall 
were studied. 

For the flow of a laminar falling film, a fully- developed flow regime 
is present where an equilibrium is established between the gravitational 
body force and the shear force at the wall. It was found that a film 
introduced at a height above or below the equilibrium height eventually 
reaches the equilibrium condition for that flow rate. In the developing 
region, the height of the free surface and the distribution of the velocity 
components change simultaneously. For plane flow along a horizontal plate, 
the film height rises up and forms a hydraulic jump where a change of the 
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flow regime from supercritical to subcritical takes place. For plane flow 
under zero gravity, the film height increases monotonically . The shape of 
the velocity profile becomes parabolic after some distance from the 
entrance but the flow never attains a fully- developed condition. 

For radial flow the height of the film and the downstream flow 
structure are determined by the flow rate, inlet height and magnitude of 
the gravitational acceleration. For larger flow rates, the film height 
decreases downstream under zero gravity whereas the film height increases 
slightly downstream when the gravitational body force is present. For 
smaller flow rates and inlet film heights, the fluid level downstream rises 
to a height above the inlet. This rise happens suddenly under normal 
gravity whereas the rise occurrs more gradually under zero gravity. The 
level of the fluid decreases slowly after the jump. This rise or jump is 
also found to be associated with a region of flow separation. The flow 
re- attaches after the separation zone and the velocity profile becomes 
parabolic downstream. 
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1.2 Introduction 


The flow of a thin liquid film is encountered in many engineering 
devices. During evaporation or condensation on a solid surface in a 
compact heat exchanger or cooling tower, spin coating in metal industries, 
and impingement of a liquid jet on a solid wall as in impingement cooling, 
a thin film is quite commonly found. Besides practical applications, the 
fluid mechanics of thin film flow is important from a theoretical point of 
view since both viscosity and free- surface effects are important in these 
flows. Moreover, the understanding of such flows under reduced or zero 
gravity is essential for proper design of heat exchangers for space 
applications. 

The falling of a thin liquid film along a plane vertical wall has been 
studied by many investigators since the turn of this century. For steady 
fully- developed laminar flow, a theoretical solution can be derived from a 
simple balance between the gravitational body force and the shear force at 
the solid wall (Bird, et al. [1]). The film height remains constant and 
the velocity profile across the film becomes parabolic in the fully 
developed region. The analysis of developing flow when a film is 
introduced at its equilibrium height is also available in the literature. 

A film falling under the influence of gravity ceases to be laminar and 
constant in thickness when the flow rate is high. Vaves tend to appear on 
the surface and the flow becomes turbulent as the flow rate is increased. 
A number of theoretical as well as experimental studies have been performed 
to understand the flow in wavy- laminar and turbulent regions. A review of 
such studies has been prepared by Faghri and Payvar [2] . This review also 
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included the experimental studies on laminar flow with constant thickness. 


A somewhat less studied problem is the spread of a liquid film over a 
plate. Vatson [3] presented results of analytical and experimental studies 
of the radial spread of a liquid jet impinging on a horizontal plane for 
laminar and turbulent flows. By using the boundary layer approximations 
for the governing equations, he derived analytical solutions using a 
similarity transformation along with the Pohlhausen integral method. His 
analysis covered the regions where the boundary layer thickness is less 
than the film height and where the film is totally engulfed by the boundary 
layer. The effects of the gravitational pressure gradient was discussed. 
The possibility of a hydraulic jump in such a flow was also anticipated. 
However, the analysis was applicable only to regions before the jump and 
the jump height could be predicted for any given location of the jump. The 
agreement between the experimental data and the analysis was satisfactory. 

Another interesting problem of thin film research is the spreading of 
the film under the action of centrifugal force as seen in a rotating 
system. An approximate analytical solution for laminar flow on a rotating 
disk was developed by Rauscher et al. [4]. An asymptotic expansion 
technique was used where the radial spread of the fluid was perturbed to 
determine the effects of convection, Coriolis acceleration, radial 
diffusion, surface curvature and surface tension. These higher order 
effects were discussed on a physical basis. The turbulent counterpart of 
the rotating disk problem was solved by Murthy [5] where the governing 
transport equations were integrated across the thickness of the film and 
the reduced equations were solved analytically. 
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The heat transfer in a thin film has also been studied by many 
investigators, particularly dealing with evaporation or condensation. The 
earliest study in this regard is credited to Nusselt [6], who developed an 
analytical solution for film condensation along a vertical isothermal wall. 
The analysis given by Nusselt has been improved by several investigators to 
include the effects of inertia, vapor pressure, etc. Sparrow and Gregg [7] 
developed a theory for rotating condensation where the removal of the 
condensate took place due to centrifugal force instead of the force due to 
gravity. Butuzov and Rifert [8] performed experiments to verify the theory 
developed in [7] . In a more recent study, Butuzov and Rifert [9] presented 
experimental as well as theoretical results for the reverse problem of film 
evaporation from a rotating disk. 

In the previous studies concerning thin liquid films, the 
investigators have tried to develop analytical models or have taken 
experimental data. Some of these models are quite approximate in nature 
and do not bring out the finer details of the flow field. A numerical 
finite- difference solution of a thin film flow is not available at the 
present time. These flows are difficult to solve by the finite- difference 
method since the geometry of the free surface is not known ahead of time, 
and the surface profile cannot be fitted in a regular Cartesian or 
cylindrical coordinate system. Moreover, none of the studies mentioned 
above has considered the flow under reduced or zero gravity, which is 
expected to be different from the flow under normal gravity. A proper 
understanding of such flows is essential in the design of space cooling 
systems. 
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The present study is undertaken to develop a general numerical 
solution procedure for thin film flows which can be applicable to both 
plane and radial systems, and to both normal and zero gravity environments. 
The numerical solution is supplemented with a one- dimensional analysis. 
The results highlight the effects of gravity for different configurations 
of the flow. 

1.3 Problem Formulation 

The equations governing the conservation of mass and momentum in a 
thin film of fluid which is newtonian with constant properties are given by 


V • V = 0 (1.1) 

-4 

9 SI = - Vp - II V 2 V + n (1.2) 

These governing equations have to be supplemented with appropriate 
boundary conditions. At the solid wall the no- slip condition exists, 

therefore, V = 0. On the free surface the shear stress vanishes which 

-4 

implies d\/dn = 0, where n is the coordinate normal to the free surface. 
Moreover, in the absence of any significant surface tension, the static 
pressure on the free surface must equal the ambient pressure. By setting p 
equal to the difference between the actual and ambient pressures, then p = 
0 on the free surface since pressure is a scalar quantity. Boundary 
conditions must also be assigned in the direction of the flow at two 
locations if the problem is elliptic. The locations usually are the 
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entrance and exit to the control volume or computation domain. At the 

entrance, the film height and velocity are prescribed: h = h^ n , and V = 

V^ n . These also determine the flow rate and the Reynolds and Froude 
numbers. At the exit the flow condition is usually unknown, so the fully 

developed condition dV/dn = 0 is prescribed. Moreover, if the film height 
or Froude number is known at the exit, this information can also be 
prescribed as a boundary condition. A boundary condition for the pressure 
has to be given at one location in the direction of the flow. In the 

numerical computation, it is convenient to prescribe this condition at the 
outlet of the flow field. The coordinate system used here is shown in Fig. 
1.1. The boundary conditions for different flow systems are listed in 
Table 1 in the component form. The flow configurations considered in the 
present investigation are shown in Fig. 1.2. They can be broadly 
classified into two groups: 

(P) Plane thin film flow 
(R) Radial thin film flow 

In the first group, three problems will be considered: 

(PI) Falling film along a vertical wall 
(P2) Plane flow along a horizontal wall 
(P3) Plane flow under zero gravity 

The first case is a classical problem where the major driving 
mechanism is the gravitational body force. This problem was used to check 
the accuracy of the present numerical scheme. The effects of introducing 
the film at a height other than the equilibrium height for a given flow 
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rate will be investigated for this case. 

For plane flow along a horizontal wall, different flow regimes may be 
present depending on the film height and fluid velocity. For plane flow 
under zero gravity, the orientation of the solid wall is immaterial. The 
flow remains the same whether the plate is vertical, horizontal or 
inclined. In this case, the flow is driven by inertial and viscous forces. 

In the second group, we will again consider two problems according to 
the presence or absence of gravity. 

(Rl) Radial flow under zero gravity 

(R2) Radial flow along a horizontal plate under normal gravity. 

As the flow spreads out radially, the cross sectional area of the flow 
increases downstream. The velocity decreases downstream due to friction 
and the change in area. The flow field changes continuously as the fluid 
moves downstream and can never become fully developed. Different phenomena 
may also happen according to the rate of the flow. 

1.4 Theoretical Analysis 

Ve consider the one- dimensional flow of a thin film where the film 

height h is a function of the radial (longitudinal, for plane flow) 

location r. Let V be the average velocity of the film in that direction 
and Q be the volume flow rate. In the case of plane flow, Q is the volume 

flow rate per unit of cross sectional area. The volume flow rate is given 

by 
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(1.3) 


Q = (2 r r) K hV 


uhPrP k = I °» for P lane s y stem 

wne e a j 1, for radial system 

Expressing the friction at the wall in terms of the friction coefficient 
gives 



It is assumed that c^ is constant in the present analysis. 

Integrating eqn. (1.2) across the thickness of the film and 
substituting the eqns. (1.3) and (1.4) results in 

3r ( 5" + gh ) = " T r (1 ' 5) 

This momentum equation has to be solved along with eqn. (1.3) to 
determine the flow field, and will be carried out in the following 
subsections. 

1.4*1 Flow Under Zero Gravity 

For a steady flow under zero gravity, the governing equations, eqns. 
(1.3) and (1.5), reduce to 

q = (2 *t) K hV = constant (1.6) 
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„ (tv c f V 2 

v a? = • r r 


(1.7) 


Eliminating h from eqns. (1.6) and (1.7) gives 



c f 


(2irr)^ dr 


Integrating this equation and substituting the inlet boundary 
condition (i.e., at r = r^, V = and h = h^) yields the solution in the 
form 



1 f r 
1 T (K + i)h 

< c f r l 

1 ' t (K nyr 



Solving for h results in the following relation 


( 1 . 8 ) 



(!-*)( 



(1.9) 


where A = 


T JK 




( 1 . 10 ) 


It can be observed in eqn. (1.9) that for plane flow where K = 0 the 
film height increases monotonically with radial location r. This increase 
in height is due to a decrease in velocity because of friction. For 
inviscid flow when no resistance is exerted by the solid wall, the film 
height remains the same at all downstream locations beginning from the 
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entrance. In an axisymmetric radial system, however, it is noticed that h 
can decrease or increase depending on the radial location. The location 
for the minimum h can be determined by differentiating eqn. (1.9) which 
gives 


( )f 0 r .in h * ( '-ir- ) 1/2 <»■“> 

( EJ )min = 2 'T 1 ^ 1 " 1 t 1 ' 12 ) 

From eqn. (1.11) we find that a minimum exists if A <: 1/2. Otherwise, the 
film height increases continuously from the entrance. 

1.4. 4 Flow on a horizontal plane under normal gravity 

This analysis is applicable for radial and planar flows under normal 
gravity with a free surface. The flow of a thin film resembles open 
channel flow, so the Froude number may be a useful parameter. Ve define 
the Froude number as 

Fr = (1.13) 

\Tgfi 

Expressing eqns. (1.3) and (1.5) in terms of the Froude number results in 
q = (2rr) K fl Fr h 3 / 2 (1.14) 

[ h ( 1 . ^ ) ] = - ^Fr 2 (1.15) 
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(a) Plane Flow 

For plane flow, K n 0, so eqns. (1.14) and (1.15) can be reduced to 
the following form: 


dFr _ 3 Fr 11 / 3 

^ 5 (Fr 2 - 1) 


(1.16) 


where 


E = 


: f g 


1/3 


(1.17) 


Integration of this equation gives 
j Fr" 8 / 3 - Fr' 2 / 3 = - K * c 


(1.18) 


To evaluate the integration constant, the Froude number must be 

specified at one location in the flow. The methodology of the present 

analysis is somewhat similar to Fanno flow analysis of one- dimensional 

compressible flow with friction. For Fanno flow, it is customary to 

evaluate the integration constant at a critical location where the Bach 

number is unity. It is assumed that a critical location exists in the flow 

* 

field. Let Fr = 1 at a location where R = R . 


Substituting this in eqn. (1.18) and rearranging gives 
R - R* = (Rr -2 / 3 - 1) - l/4(Fr~ 8 / 3 - 1) (1.19) 
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Equation (1.19) is a double- valued function as shown in Fig. 1.3. The 
two stems of the function represent subcritical and supercritical flows 
where the Froude number is less and greater than unity, respectively. 


Since two solutions exist at any location, the possibility of a sudden 
jump from supercritical to subcritical flow exists. The opposite is not 
true since that would violate the second law of thermodynamics. The height 
of the film before and after the jump can be related by the conservation of 
mass and momentum across the jump. For plane flow this is given by 

H7 = j [ J 1 - 8 Frj 2 - 1 ] (1.20) 

where subscript 1 indicates conditions before the jump and subscript 2 
indicates conditions after the jump. 


(b) Axisvmmetric Radial Flow 

For radial flow, K = 1, so eqns. (1.14) and (1.15) can be transformed 
into the following form 


dFr _ Fr (2 + 
2E (Fr 2 


Fr 2 ) 3Fr n / 2 R 2 / 3 
- 1) 2 (Fr 2 - 1) 


( 1 . 21 ) 


where 


R = 



( 1 . 22 ) 
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The radial distance has been non-dimensionalized in terms of other flow 
parameters. An analytical solution for eqn. (1.21) is not possible, so a 
numerical integration was performed using the Euler method. It can be 
noticed that the equation is singular at Fr = 1. Therefore, the critical 
condition cannot be directly applied as a boundary condition in the 
numerical solution. To avoid this singularity, the equation can be 
expanded around the singular point and the solution can be found at a short 
distance from the singular point from the lowest- order expansion. The 
numerical integration was then carried out beginning from a short distance 
away from the singular point, where the solution is already known. The 
solution is shown in Fig. 1.4. 

* 

It should be noted that the critical radius, R , appears as a 
parameter. Also, the shape of the curve in somewhat different from that 
for plane flow. This is because the flow spreads out as it propagates 
radially downstream. The double- valued nature of the solution still 
remains in radial flow which indicates the possibility of a hydraulic jump. 

1.4.3 Characteristic Behavior of the Flow 

Since the equations of transport in a thin film are somewhat similar 
to those for one- dimensional compressible flow, it may be useful to analyze 
the characteristic behavior of the flow. 


The conservation equation in its time- dependent form can be written as 


5U 1 d 
M + pr 3r 


(r«E) 


= H 
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where 



v, 


hV 

u = 

n 

hV 

, E = 

h(V 2 + 


and H = 


o 

JP 


where is the surface shear stress. These are two first- order ordinary 
differential equations in t and r with two dependent variables, h and V. 

Expressing the equations in the characteristic form results in the 
following equation 


*t + C R x = S 


where 


v * 'CPI 

v - c-pIk; 


= Rieraann invariants 


C = 


v + y ghj_ 

V - V ghj 


= Wave speed 


It can be seen that the first invariant always propagates downstream. The 
second invariant, however, propagates downstream for supercritical flow (Fr 
= V/V ghj > 1) and propagates upstream for subcritical flow (Fr < 1). This 
implies that both V and h must be prescribed upstream for supercritical 
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flow while only one is prescribed upstream and one downstream for 
subcritical flow. This analysis will assist in determining the appropriate 
location of the boundary conditions for the computations addressed in the 
following sections. 

1.5 Numerical Solution Procedure 

The governing transport equations along with the appropriate boundary 
conditions were solved numerically using a finite- difference scheme. Since 
the free surface geometry cannot be handled very well with a regular 
rectangular or cylindrical coordinate system, a boundary- fitted curvilinear 
coordinate system had to be used. In this system, the free surface of the 
film was used as one of the boundaries of the control volume. 

A curvilinear system can be either orthogonal or non- orthogonal 
depending on whether the faces of the control cells are orthogonal to each 
other or not. The orthogonal system has the advantage of simplicity 
compared to the non- orthogonal system. In either system, the vectorial 
form of the governing equations [i.e., eqns. (1.1) and (1.2)] can be 
written in terms of components and can be discretized to determine the 
finite- difference equations. In most of the computations presented, the 
coordinate system was non- orthogonal. 

Vithin the range of the general non- orthogonal coordinate system there 
exist several options in formulating the equations. These options arise 
from the freedom available in the choice of velocity components and their 
direction with reference to the coordinates. Thus, velocity and force 
vectors can be resolved either into their Cartesian, covariant or 
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contravariant components. Moreover, the problem can be solved in a 
physical domain or transformed into a domain where the grid cells are 
rectangular and other physical quantities are non-dimensional or reduced in 
dimension. Although all these options are obviously equivalent to each 
other from the physical point of view, they are substantially different as 
far as numerical treatment is concerned, each presenting its own problems. 
Some of these problems are highlighted in the paper by Galea and Markatos 
[10] and are discussed briefly below. 

The resolution of the velocity and other vectors into Cartesian 
components is possibly the simplest approach, but this leads to an 
inaccurate numerical formulation, in particular for the pressure gradient 
term driving the flow. The use of contravariant components makes the 
formulation of convective terms straightforward, but still does not solve 
the problem of accurately formulating the pressure gradient and does not 
help with respect to properly assigning the value of the velocity convected 
into the distorted grid cell. The use of covariant components solves the 
two problems above, but leads to more complicated calculations for the 
convective terms. In the present study, the problem was solved in the 
physical domain where covariant velocity components were used. 

The grid system used can be considered as a distorted version of the 
usual orthogonal Cartesian grid system in which grid lines and control 
cells are stretched, bent and twisted in a arbitrary manner, subject to the 
cells retaining their topologically Cartesian character. This means that 
grid cells always have six sides and eight corners in a three-dimensional 
case. 
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As shown in Fig. 1.1, the local coordinates are defined along lines 
joining adjacent cell centers. The z-axis was taken in the stream line 
direction and y-axis in the direction across the film. The resolutes of 
the velocity vectors in the y and z directions are v and w, respectively, 
and can be defined as 


v 



j 


w = 


V 


• k 


Here j and k are unit vectors in the direction of the coordinate axes. In 
general, the resolutes are not the same as the velocity components in these 
directions, but can be related to them by geometrical factors. 

The finite difference equations were derived by the application of the 
conservation principle of mass and momentum to the grid cells. The 
transport processes for each cell are convection and diffusion. Moreover, 
there may be a momentum or mass source within the cell. The mass flux 
across a cell boundary was computed exactly from the scalar product of the 
velocity vector and the vector representing the area of the cell face. 
Note that this can be written out in terms of velocity resolutes and 
geometrical factors including angles between cell faces. In the 
calculation of convection across a cell face, special attention was given 
to the change of orientation of the coordinate axes from cell to cell and 
the curvature of a cell face. These resulted in extra terms in the 
calculation of convection. However, the representation of convection was 
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exact and did not involve any approximation due to the non- orthogonality of 
the grids. 

The calculation of diffusion is somewhat more complicated than 
convection. The diffusion flux was calculated assuming the coordinate 
system to be locally orthogonal. This obviously neglects cell curvature 
and non- orthogonal orientation and may incorporate a substantial amount of 
error where the process is primarily diffusive. However, in the thin film 
calculation this approximation should not introduce severe inaccuracies, 
particularly when the film enters the control volume with a reasonably high 
velocity. 

The relative importance of convection and diffusion at each cell was 
determined from the magnitude of the local Peclet number. A hybrid 
difference scheme demonstrated by Patankar [11] was used. The calculation 
of the momentum source due to the pressure gradient and that due to the 
gravitational body force could be accomplished without any approximation 
for non- orthogonality . 

The grid generation was achieved in two steps. First, the grid cells 
were formed by algebraic interpolation between the boundary points. This 
provided an approximately equal volume for each control cell. The 
boundaries for the interior cells were then smoothed to make the cell faces 
more orthogonal to each other. This operation resulted in a better 
representation of diffusion in the flow field and more accurate 
computations. The details of the formulation in a body- fitted coordinate 
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system and the generation of grid cells are described in the work by 
Spalding et al. [12]. 

The number of grids in the y-z plane was determined from a series of 
computations with different grid sizes. A sample of this investigation is 
shown in Fig. 1.5 where the free surface profile for radial flow under zero 
gravity is plotted. It can be seen that the free surface profile is 
somewhat sensitive to the grid size. The solution obtained by using 20 x 
20 grids is somewhat different from that obtained with 40 x 20 grids in z-y 
directions. The differences become smaller when the 40 x 20 solution is 
compared with the 50 x 25 solution. Finally, the difference completely 
disappears when the solution corresponding to 50 x 25 grids is compared 
with the 55 x 27 solution. Therefore, 50 grids in the direction of flow 
and 25 in the direction across the film was found to be adequate and was 
used for all radial flow computations. For plane flow, it was found that 
the solution corresponding to 40 x 20 grids in the z-y plane precisely 
predicts the friction coefficient and velocity profile in a falling film 
system. Therefore, further refinement of the grids appeared to be 
unnecessary and all computations for plane flow were carried out using 40 x 
20 grids except for case P2, which involved a complex free surface profile. 
For this case 50 x 25 grids were used. 

The flow field was solved by using SIMPLEST algorithm as discussed by 
Spalding [13] . One special feature of this algorithm is that in the 
discretized form of the momentum equation, the convection terms are lumped 
together with the source term. This results in a faster convergence for 
some flow conditions. The algorithm works in an iterative manner where the 
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continuity equation is transformed and used as a pressure correction 
equation. The computation starts by guessing a pressure field. This is 
used to determine velocity components from their corresponding momentum 
equations. The modified continuity equation is then used to determine the 
amount of pressure correction. The guessed pressure, the amount of 
pressure correction and the solution corresponding to the momentum 
equations are then assembled together to give the flow rate and pressure 
field for that step. The new pressure serves as a guess for the next step. 
The solution proceeds across the control volume until the normalized 

- fi 

residual for each equation was approximately 10". The above 

finite- difference formulation and solution procedure is incorporated in 
computer program PHOENICS that was used in the present study. In the free 
surface flow discussed here, both the zero- shear condition and the p = 0 
condition at the free surface need to be satisfied. These two conditions 
cannot be simultaneously given at a boundary. On the other hand, the free 
surface geometry, which is unknown in the problem has to be given before 
solving the flow field by a finite- difference method. To avoid this 
difficulty, an iteration scheme has been adapted as described below. 

(1) Guess a free surface profile. 

(2) Solve the flow field completely for that profile using the 
zero- shear condition on the free surface boundary. 

(c) Find the pressure distribution on the free surface and calculate 
its deviation from an ideal zero- pressure free surface. The 
measure used here is the RMS (root- mean- square) error. 

(d) Calculate and reduce the RMS error on the free surface by 
successive alteration of the surface profile. 
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(e) The results with the minimum error give the required final 
solution. 

The accuracy of the computation depends somewhat on the assumption of 
the form of free surface. In the results presented here, hyperbolas with 
two or more degrees of freedom were used to represent the computation 
domain. 


1.6 Results and Discussion 
1 . 64 Plane Film 

The flow of a plane film along a vertical wall under the influence of 
gravity is a classical fluid mechanics problem where an analytical solution 
is available for fully- developed laminar flow. In the fully- developed 
region, the film height remains constant and the velocity profile has a 
parabolic appearance. Numerical computations using the present methodology 
were performed for a film which has already reached the fully- developed 
condition. A parabolic velocity profile with the same shape as given by 
the analytical solution was used for the incoming fluid. The Reynolds 
number for the film was Re = 50 where water was used as the fluid. It was 
found that for the entire domain, the velocity profile across the film 
remains about the same and the friction coefficient was equivalent to that 
of the analytical solution. In the present investigation, the friction 
coefficient is defined as 


c f = 




(1.23) 


22 



where ¥ is the average velocity of the fluid across the film. This 
definition is applicable for both plane and radial flow. For plane flow 
with constant thickness, V is a constant because of continuity. 


The distribution of pressure on the free surface was also computed and 
the RMS error was determined to be 0.18 Pa. This non- zero figure of RMS 
error may be attributed to the inaccuracy associated with the 
representation of flow field by finite- difference equations. This amount 
of error cannot be avoided in a numerical solution, so it provides a lower 
bound on the error that can be expected in thin film computations. 

The developing flow of a falling film when introduced at a height 
equal to, or above or below the equilibrium height was also investigated. 
The flow conditions are summarized in Table 2 and the results are shown in 
Table 3 and Figs. 1.6- 1.9. When the film enters the control volume at the 
equilibrium height, the height remains the same and the development of the 
velocity profile from uniform to parabolic occurs as the flow moves 
downstream. When the film enters with a height other than the equilibrium 
height, a gradual adjustment of height takes place until the flow reaches 
the equilibrium height. The adjustment of the free surface and the 
development of the velocity profile occur simultaneously in this flow. To 
model the free surface, we have assumed a profile of the form 


h = 



( 1 + I ) an 
6 


for z < Zj 
for z > Zj 


(1.24) 


The first part of this profile provides the variation of the free surface 
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height in the developing flow region and the second part gives the height 
after the adjustment is completed. The downstream location where the free 
surface adjustment is complete is z^. 

For the computation with h^ n = 8, we find that the RMS error on the 
free surface condition is 1.2 Pa. The uniform entrance condition provides 
a higher RMS error than the parabolic entrance since a developing flow 
region is present where the pressure has to conform with the flow 
development. An error of this order may be acceptable for a free surface 
computation since no variation of height is expected to take place when a 
film is introduced at the equilibrium height. The distance required for 
flow development was found to be about 5 times the film thickness in this 
case. A flow is defined to be fully- developed when the wall shear stress 
is within 2 7. of the final equilibrium value. This definition is similar 
to that given by Kays and Crawford [14] for developing flow in closed 
conduits. 

As evident from Fig. 1.6 the adjustment of film height takes places 
for a length of 0.006 m (i.e. , approximately 105) when a film is 
introduced at a height 207. more than the equilibrium height. In this 
situation the velocity profile also becomes fully -developed at the same 
location. A film introduced at a height 207. lower than the equilibrium is 
found to take a shorter distance for the adjustment of the free surface and 
the velocity profile. The distribution of the free surface pressure is 
shown in Fig. 1.7. It can be noticed that most of the pressure variation 
takes place in the developing flow region and its value becomes very close 
to zero when the flow becomes fully developed. The RMS error in pressure 
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is small in all cases which indicates that modeling the free surface with 
an equation of the form of eqn. (1.24) is reasonably accurate. 

The variation of the friction coefficient along the length of the film 
is shown in Fig. 1.8. The distributions for h^ n = 5 and h^ n = 0.85 are 
very close to each other, but for h^ n = 1.25 the variation is significantly 
different. The friction coefficient first decreases to a value lower than 
that for the equilibrium condition and then rises up to the equilibrium 
value. This is due to two counter- acting phenomena that affect the film in 
the developing flow region in this case. The wall friction propagates 
outward as the boundary layer develops beginning from the entrance point 
which tends to reduce the velocity of the fluid. The thickness of the 
film, however, decreases and tends to increase the fluid velocity due to 
the area available for the flow. Since the first effect starts from the 
wall, it is more dominant in the earlier part of the flow development and 
then the second effect takes over in the region downstream. It can also 
mentioned that when the film is introduced below the equilibrium height, 
the increase in the film height and the propagation of the shear tend to 
reduce the film velocity, so the behavior is not analogous to the case when 
the film is introduced above the equilibrium height. In all situations, a 
plane falling film eventually attains a fully- developed flow. This is 
confirmed by comparison of the velocity profile in a location near the exit 
as shown in Fig. 1.9. The variation of the velocity and the friction 
coefficient are identical in all situations. 

The flow of a plane film was also investigated for a horizontal 
orientation of the plate where the gravity acts across the thickness of the 
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film instead of in the direction of the main flow. The flow conditions 
used here are listed in Table 2 and the corresponding free surface geometry 
is described in Table 3. Figs. 1.10-1.13 present the computational 
results. 

In Fig. 1.10 it can be seen that after the film enters the control 
volume, the liquid film height increases rapidly within a short distance 
until it is about 4 times the inlet height. The film height then levels 
off slowly as the flow moves downstream. Close to the exit the height 
decreases to accommodate the outlet boundary condition. In this case, the 
inlet Froude number is 3.0 which means that the flow enters the control 
volume at a supercritical condition. The transition of the flow from 
supercritical to subcritical takes place with a rapid change in fluid 
level. Once the subcritical condition is established, the change in height 
is basically due to the deceleration of the flow by friction at the wall. 
Figure 1.10 also shows the results obtained by Thomas [15] for a 
one- dimensional model of the flow under the same flow conditions. This 
result was obtained by numerical integration of the subcritical and 
supercritical streams proceeding from the end points and matching the 
conditions at the jump location. The one- dimensional solution seems to be 
quite different from the two-dimensional flow considered here, particularly 
in the subcritical region. It appears to be quite strongly influenced by 
the downstream boundary condition whereas in the two-dimensional case the 
influence is less dramatic. A computation using the present numerical 
scheme where the free surface was forced to take the form obtained by 
Thomas [15] showed a RMS error of 21.62 Pa on the free surface pressure, 
which is about 10 times the RMS error in our final solution, as listed in 
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Table 3. Therefore, the free surface profile predicted by the 2-D solution 
is likely to be much more accurate than that predicted by the 1-D solution. 

The distribution of the flow velocity across the thickness of the film 
is shown in Fig. 1.11 for two locations on the plate. The velocity profile 
is parabolic in nature. Even though the average velocity decreases 
downstream, the local velocity may increase or decrease depending on the 
location across the film thickness. In the numerical results, it was also 
seen that there is a separated flow region around the location of the jump. 
To cope with the rapid rise of the fluid level, flow separation occurs near 
the wall, which reattaches again at a downstream location. The profiles 
shown here correspond to locations beyond the reattachment point. Inside 
the recirculating flow region, the fluid moves in a direction opposite to 
the main flow at locations adjacent to the wall. Even though the velocity 
is positive at all locations after reattachment, the acceleration of the 
fluid near the wall continues for a distance as seen in Fig. 1.11. 

The shear stress exerted by the wall and the corresponding friction 
coefficient are shown in Figs. 1.12 and 1.13, respectively. The wall shear 
stress decreases rapidly after the entrance, which is due to the 
simultaneous adjustment of the film height and the development of the 
velocity profile from uniform entrance conditions to a two-dimensional 
form. Both effects aid to a decrease in the shear stress. The shear 
stress is negative in the recirculating flow region and becomes positive 
again after reattachment. The two locations corresponding to zero shear 
are the separation and reattachment locations. In Fig. 1.12, it can be 
also noted that the wall shear becomes constant at an intermediate location 
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of the flow and rises again near the exit. The increase in the wall shear 
stress and the friction coefficient are due to an increase in flow velocity 
near the wall as the exit is approached. 


The remaining part of the investigation concerning plane film flow is 
the flow of a film under zero gravity. In this situation the orientation 
of the plate become immaterial and an identical flow condition is achieved 
if the plate is horizontal, vertical or inclined. Consider a situation 
where the film is introduced at a height equal to the equilibrium 
fully- developed flow in a fallings film system. In the absence of gravity, 
the flow is acted on only by viscous and inertial forces and the film 
height is expected to increase downstream. To model the free surface, we 
assume a profile of the form given in Table 2. The flow conditions are 
also listed in the same table. The results are shown in Table 3 and Figs. 
1.14-1.18. 


Figure 1.14 shows the variation of the film height with distance which 
increases monotonically. The figure also shows the analytical solution 
derived in the previous section. The analytical solution requires the 
specification of a friction coefficient. In the present investigation it 
was taken from the numerical solution instead of assuming it to be constant 
throughout the region. The comparison between the analytical and numerical 
solutions appears to be good in most regions of the flow. 

The variation of the free surface pressure is shown in Fig. 1.15. 
Analogous to the falling film flow, most of the pressure variations occur 
in regions close to the entrance where the development of the velocity from 
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a uniform to parabolic profile takes place, which is analogous to the 
falling film flow. However, in contrast to the falling film, the velocity 
keeps changing as the flow moves downstream and never attains a 
fully- developed situation. To illustrate this point, the velocity profile 
at three different locations on the plate are plotted in Fig. 1.16. Even 
though the profile becomes parabolic in nature, the magnitude of the 
velocity decreases continuously. The shear stress exerted by the solid 
wall is shown in Fig. 1.17 and the corresponding friction coefficient is 
plotted in Fig. 1.18. It can be noticed that the shear stress decreases 
continuously as the flow moves downstream, whereas the friction coefficient 
has a minimum at an intermediate location and then increases downstream. 
The largest variations of the shear stress occur close to the entrance due 
to the development of the velocity from a uniform to a parabolic profile. 
After the velocity profile is parabolic, the slight reduction in the shear 
stress is due to the deceleration of the flow, which is small compared to 
the reduction of the average velocity in regions away from the entrance. 
This results in an increase in the friction coefficient. 

1.6.2 Radial Film Flow 

The system where a fluid is introduced at the center of a circular 
horizontal plate and spreads uniformly in all radial directions was also 
studied. Two combinations of flow rates and initial film heights were 
chosen and water was used as the fluid. These cases were studied for 
zero-gravity (g = 0) and normal gravity (g = 9.81 ra/s ) situations. The 
effect of critical outflow as may happen at the exit of the plate was 
studied for one case. The flow parameters for different cases studied here 
are shown in Table 4. The surface profile and corresponding RMS error for 
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the surface pressure for different cases are listed in Table 5. The 
constants appearing in each profile equation were optimized to give the 
minimum possible error in the free surface pressure. The values appearing 
in Table 5 are optimum values corresponding to each flow configuration. It 
can be noticed that all RMS errors are on the order of 1 Pa, so the 
solutions may be considered to be acceptable. 

Figure 1.19 shows the variation of the film height with the radial 

distance for cases RlA and R2A. These cases correspond to the same flow 

rate and inlet film height but the gravitational acceleration for RlA is g 

o 

= 0.0 and for R2A is g = 9.81 m/s . The one- dimensional analytical 

solutions, as discussed in a previous section, corresponding to these cases 
are also shown in the figure. The analytical solution requires the 
specification of a friction coefficient which was taken to be the average 
friction factor in the numerical solution. The comparison between the 1-D 
and 2-D solutions appears to be quite reasonable. In the derivation of the 
one- dimensional solution it was assumed that the film has a uniform 

velocity across its thickness, whereas in reality the film has a parabolic 
velocity profile as seen in Fig. 1.20. This assumption, along with the 
assumption of a constant friction coefficient may be the primary reasons 
for the difference between the 1-D and 2-D solutions. 

Comparing the film height for the case RlA with that for the case R2A, 
it can be seen that gravity has very significant effect in determining the 
flow behavior of the film at this flow rate. In the absence of gravity, 
the film height gradually decreases in the downstream direction, whereas 
when g = 9.81 m/s , the film height slightly increases as the flow proceeds 
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downstream. In the absence of gravity, there is no static pressure 
gradient across the thickness of the film and the momentum balance in the 
film is solely due to inertial and viscous effects.. The pressure variation 
is significant for a film of moderate thickness which tends to flatten out 
the free surface. 

The variation of velocity across the film thickness is shown in Fig. 
1.20 for two different locations. The profile remains parabolic but its 
magnitude decreases downstream due to the increase in area available for 
the flow. The decrease is less dramatic for the zero gravity case since 
the film height also decreases downstream. The distribution of free 
surface pressure for the two cases are plotted in Fig. 1.21. Figures 1.22 
and 1.23 show the variation of the wall shear stress and the friction 
coefficient at this rate of flow, respectively. The largest variations of 
these quantities appear to be in regions close to the inlet. This is 
primarily an outcome of the development of the velocity profile from 
uniform at entrance to parabolic in downstream locations. For case E2A, 
the wall shear decreases downstream, which is expected since the velocity 
of the fluid decreases and consequently encounters less resistance from the 
wall. In case R1A, however, the wall shear stress increases slightly at 
locations far downstream. This is due to the local acceleration of the 
flow in the region next to the wall. It may be pointed out that even 
though the velocity profile is parabolic, the actual shape of the parabola 
changes along the plate. This is caused by the increase of the local 
velocity next to the wall. The increase of the friction coefficient for 
cases RlA and R2A is primarily caused by a decrease of the average velocity 
as the fluid moves downstream. 
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Figures 1.24-1.28 demonstrate the behavior of the flow when the film 

is introduced at a smaller height and flow rate. As in the previous two 

cases, case E1B and R2B correspond to identical flow and boundary 

conditions, but the gravitational acceleration in case RIB is g = 0.0 and 

o 

in case R2B is g = 9.81 m/s . Case R2C is an extension of case R2B where 
the film at the exit is forced to have a Froude number of Fr = 1. When the 
radius of the plate is finite and ends at the exit of the computation 
domain, the fluid spreads over the plate and then free falls over the edge 
of the plate. For a free- falling exit situation, a subcritical flow (Fr < 
1) has to transform to a supercritical flow (Fr > 1) at the discharge 

location. To incorporate this effect, the Froude number at the edge of the 
plate is set to Fr = 1, which corresponds to critical flow. Cases RIB and 
R2B does not take into account this effect and the plate is assumed to 
extend somewhat beyond the end of the computation domain. 

The film height for the three different cases are shown in Fig. 1.24. 
The analytical solution corresponding to case R2B is shown here for 

comparison. At this condition, the fluid rises up and forms a hydraulic 
jump close to the entrance. The shape of the free surface and its maximum 
height depend on the gravitational body force and exit boundary condition. 
In case RIB, which corresponds to a situation where g = 0.0, the increase 
of the fluid level is about 107# more than the other two cases when g = 9.81 
m/s . The presence of gravity, which causes a pressure gradient across the 
thickness of the film, tends to flatten out the free surface profile. It 

is also evident that a sudden rise of the flow occurs when gravity is 

present which contrasts the gradual increase of the fluid level when g = 
0.0. After the jump, the height of the film slowly decreases as the fluid 
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proceeds further downstream. In case E2C, when the outlet is constrained 
to have Froude number of Fr = 1, the film maintains about the same 
thickness for the most regions of the plate and suddenly drops very close 
to the exit. In Fig. 1.24, it can be also noticed that the height and 
location of the jump is the same for both cases R2B and R2C. This may 
suggest that the jump height is determined by the flow rate and Froude 
number before the jump, rather than conditions at the outlet. The 
conditions of the flow after the jump, however, depend on the exit 
conditions as evident from the shape of the free surface. This observation 
is consistent with an open channel flow where the flow rate and Froude 
number before the jump determine the height after the jump, and the fluid 
motion downstream is subcritical and affected by the outlet condition. 

Figure 1.24 also shows the analytical solution corresponding to 
one- dimensional radial flow in the presence of gravity. The solution 
requires the specification of the friction coefficient and one boundary 
condition. The friction coefficient used here is 2.0, which is the average 
value after the jump as obtained in the numerical solution corresponding to 
case R2B or R2C. When the inlet condition is used as the boundary 
condition, the analytical solution predicts a jump to a height close to 
that of the numerical solution and the film height diminishes slowly 
thereafter. The figure also shows the film height when the exit boundary 
condition corresponding to case R2B is used. The solution shows that the 
film height gradually increases in the upstream direction, but no drop in 
fluid level is obtained near the inlet. This behavior of the analytical 
solution reaffirms the previous observation that the exit boundary 
condition in the subcritical region cannot influence the supercritical flow 
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before the jump. The numerical solution is between the two analytical 
solutions. The one- dimensional solution does not account for the details 
of the behavior of the flow like separation after the jump and the 
parabolic velocity profile thereafter. The difference between the 
analytical and numerical solutions may be attributed to these fine details 
of the flow field. 

The velocity profile across the film at two radial locations are shown 

in Fig. 1.25. The velocity profile appears to be parabolic at each 

location both in the presence and absence of gravity. The magnitude of the 

velocity at z/z m = 0.8 is less than that at z/z ffl = 0.4 due to the larger 

area available for the flow. At a particular location, the velocity is 

o 

greater for the case of g = 9.81 m/s than the corresponding zero-gravity 
case. This is because of the smaller film thickness in the presence of the 
gravitational field at this flow rate. The free surface pressure for the 
three cases studied here are shown in Fig. 1.26. The pressure at 
intermediate locations become very close to zero. There are some 
fluctuations close to the inlet and the exit for the case R2C. This is due 
to the rapid change of the flow in these regions to account for the imposed 
boundary conditions. The assumed form of the free surface profile may not 
be very precise in these regions. A more complex profile assumption with a 
larger number of degrees of freedom may smooth out these pressure 
fluctuations. However, the profiles used in the present investigation are 
reasonably simple to handle in computation and capture the most significant 
structures of the flow and yield solutions that are physically realistic 
and accurate enough for engineering applications. Therefore, it is 
preferable to use a rather simple free surface profile rather than making 
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it very general to fit all regions of the flow. It should be also 
mentioned that the present free surface profiles were also arrived at after 
a number of trials with different parabolas, hyperbolas and polynomials for 
each case studied and the form selected was found to be the most accurate 
and general. 

The wall shear stresses and friction coefficients for cases RIB, R2B 
and R2C are plotted in Figs. 1.27 and 1.28. In case RIB, which corresponds 
to g = 0.0, the shear stress is positive and high near the entrance. At a 
location downstream, it becomes negative and then rises up to positive 
values again and diminishes slowly thereafter. The negative value of shear 
stress is due to separation of the flow after the jump. It becomes 
positive again when the flow reattaches, and assumes a parabolic profile. 
The small decrease in the magnitude of shear stress as the flow moves 
further downstream is due to the gradual expansion of the flow because of 
the area change. In cases R2B and R2C, the shear stress starts with a 
negative value, rises up to a positive value at a location downstream and 
slowly diminishes thereafter. In these cases, the jump occurs immediately 
downstream from the entrance and results in flow separation in that region. 
It may be noticed that the value of shear stress is slightly different near 
the entrance for cases R2B and R2C. In case R2B a uniform velocity profile 
was assumed at the entrance, whereas in case R2C a parabolic velocity 
profile was assumed. Even though both profiles become distorted due to 
flow seperation, the magnitude of the wall shear stress is somewhat 
affected by the entrance condition. The separation of the flow results in 
a very complex distribution of the friction coefficient near the entrance 
as seen in Fig. 1.28. After the flow reattaches to the surface, the 


35 


friction coefficient gradually increases downstream. It may be recalled 
that the friction coefficient in the present investigation is defined in 
terms of the local velocity. The gradual decrease of the local average 
velocity as the fluid spreads downstream overpowers the loss of wall shear 
stress and increases the friction coefficient downstream. 

1.7 Conclusions 

The present study developed a numerical solution procedure for the 
computation of plane or radial thin film flow in a normal or zero gravity 
environment. An analytical solution was also derived for a one- dimensional 
approximation of the flow. A reasonable agreement between the numerical 
and analytical solutions was obtained for most flow configurations 
considered here. 

Five different flow systems were studied: (a) plane falling film, (b) 
plane flow on a horizontal plate, (c) plane flow under zero gravity, (d) 
radial flow on a horizontal plate, and (e) radial flow under zero gravity. 
For plane falling films, an equilibrium fully- developed region is present 
where the gravitational body force is balanced by the viscous shear force 
at the wall. The film, whether introduced at the equilibrium height or at 
a height above or below the equilibrium, eventually comes to this 
fully- developed condition. In the fully developed region the numerical 
velocity profile, which is parabolic, matched exactly with the analytical 
solution irrespective of the entrance condition of the film. When a film 
is introduced at its equilibrium height and with the same velocity profile 
as that of fully- developed flow, no further change in the velocity 
distribution is found to happen as the film proceeds downstream. This 
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gives confidence in the stability and accuracy of the present numerical 
scheme. The error on the estimation of the free surface pressure in this 
case is the minimum that cannot be avoided in a numerical computation. An 
estimate on the error bound of the free surface pressure was developed from 
the calculation of developing flow and it was concluded that an error of 
the order of 1 Pa may be acceptable. The length of the developing region 
was found to be small. Both the free surface height and velocity profile 
appeared to arrive at the equilibrium condition within 10 times the 
equilibrium film thickness. The length was found to be relatively larger 
when the film begins with a height^ above the equilibrium height. 

For plane flow on a horizontal plate it was found that the film height 
increases rapidly and forms a hydraulic jump. The flow regime changes from 
the supercritical entrance condition to a subcritical downstream condition. 
The analysis of the flow indicated that both supercritical and subcritical 
flow proceed towards a critical flow when friction is present. The sudden 
rise of the fluid level resulted in flow separation adjacent to the solid 
wall. The flow reattached again at a downstream location and developed a 
parabolic velocity profile. The location of separation and reattachment 
could be determined from a plot of shear stress at the wall, which vanishes 
at these locations. 

For a plane flow under zero gravity it was found that the film 
thickness monotonically increases as the flow moves downstream. The 
velocity profile is parabolic except for regions very close to the 
entrance. The shear stress at the wall decreases as the flow moves 
downstream. The friction coefficient was computed in terms of the local 


37 


average velocity and was found to increase after coining to a minimum at an 
intermediate location on the plate. 

For radial flow, computations were performed for two sets of flow 

rates and inlet heights. The behavior of the flow was found to depend 

significantly on the flow rate and magnitude of the gravitational 

acceleration. At the higher flow rate, it was found that the film 

decreases in thickness downstream under zero gravity whereas the thickness 

o 

increases downstream when the gravitational acceleration is g = 9.81 m/s . 
The velocity profile was parabolic in most regions, even though the shape 
of the profile changed as the flow moved downstream. The change here was 
more severe than plane flow since the area available for the fluid 
increased downstream. For zero gravity it was found that the wall shear at 
locations far downstream increased due to the local acceleration of fluid 
near the wall. 

When the film thickness at the inlet is low and the Froude number is 
greater than unity, the flow jumped up within a short distance from the 
entrance and remained subcritical downstream. The jump happened at the 
same location and height whether or not the Froude number at the exit was 

Fr = 1. The flow after the jump was affected by the outlet condition. 

For the same flow condition when the gravitational acceleration is reduced 

to zero, it was found that the rise in fluid level still occurs, even 

though the rise is more gradual and to a greater height than the 
corresponding jump case. The jump or rise in fluid level was also found to 
be associated with flow separation. The flow was found to reattach to the 
surface and form a parabolic profile after the jump. After the 
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reattachraent, the wall shear stress slightly decreased as the flow moved 
downstream, irrespective of whether or not gravity was present. For radial 
flow, the friction coefficient was found to increase downstream in regions 
somewhat away from the entrance and flow separation., 
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Nomenclature 

A Optimization parameter in free surface equation, flow parameter 
(defined by eqn. (1.10)) 

an Optimization parameter in free surface equation 

bn Optimization parameter in free surface equation 

C Optimization parameter in free surface equation, wave speed 

c^ Friction coefficient 

Fr Froude number, V/V gh 

g Acceleration due to gravity 

h Film height 

h^ n Inlet film height 

h Qut Film height at exit 

p Static pressure 

t Time 

K Exponent for eqn. (1.3) 

Q Flow rate 

r Radial coordinate 

* 

r Radial location for critical flow 
Re Reynolds number, 4Q/i/ 

R Riemann invariant, dimensionless radial coordinate (eqns. (1.17) and 

( 1 . 22 )) 

* 

R Critical radius (dimensionless) 

v Normal component of velocity 

V Velocity vector 

w Component of velocity along the plate 

V Average velocity along the plate 
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y Coordinate normal to main flow direction 

ym Maximum value of y 

z Coordinate along the plate 

zm Maximum value of z 

Zj Optimization parameter in free surface equation 
Optimization parameter in free surface equation 

Greek symbols 

p Density 

p, Dynamic viscosity 

r Shear stress at wall 

w - 

6 Equilibrium film thickness 

v Kinematic viscosity 
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Case PI: 


Cases P2 
and P3: 


Cases HI 
and R2: 


Table 1 : Boundary Conditions for Plane and Radial Flow 


at y = 0: 
at y = h: 
at z = 0: 

at z = L: 


v = 
dv 


w = 0 

= 0, v = 0, p = 0 

w = V. , for uniform entrance 
in’ 

w = 1,5 V in [ 2 (f) ' (f) 2 ]’ for Parabolic 

entrance 


p = 0 


at y = 0 
at y = h 
at z = 0 
at z = L 


v = v = 0 

|^ = 0, __v = 0, p = 0 

w = V. 
m 

P = Pg (h-y) 


at y = 0: 
at y = h: 
at z = r in : 

at z = r out : 


v = w = 0 


^-0, v-0, p-0 

w = V^ n , for uniform entrance 

w = 1.5 V in [ 2(1) - (f) 2 ], for 


parabolic 

entrance 


P = Pg (h-y) 
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General data: 


Case P1A: 


Case P1B: 


Case PIC: 


Case P2: 


Table 2: Flow Conditions: Plane Flow 

p = 913.026 kg/m 3 
v = 7.432 x 10' 6 m 2 /s 


h. = 0.000595 m 
in 

V in = 0.1561 m/s, g = 9.81 m/s 2 

Free surface profile: h = 6, where 6 = 0.000595 m 

h in = 0.000714 m (1.2 S) 

V in = 0.1301 m/s, £ = 9.81 m/s 2 

Free surface profile: f b b in ^ + j) For z < Z 1 

*• h = 6, for z > z^ 

h in = 0.000476 m (0.8 6) 

V in = 0.1952 m/s, g = 9.81 m/s 2 

Free surface profile: f b b in ^ + For z < Z 1 

*■ h = 6, for z > Zj 


h in = 0.0009674 m, Fr 4n = 3.0 


m 


V in = 0.2922 m/s , z out = 0.1445 m 

Fr out = 1 -°» h out = °- 002012 g = 9.81 m/s' 


Free surface profile: 


h = h in + l)» for z < z i 
h = Ch in [2 -(1 + |) bn ], for Zj<z<z 2 
( z - z o) 

h = M. = 2 


a 2 + (z 


out 


z o) 


(h 


out 


b z = z 9 \ for z > z f 
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Case P3: h^ n = 0.000595 m 

V in = 0.1561 m/s, g = 0.0 

‘Z an 

Free surface profile: h = h^ n (1 + -j) 
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Table 3: Summary of Results: Plane Flow 


RMS Error in 

Problem Optimum Profile Pressure (Pa^ 


Case P1A (Parabolic Inlet) 

- 

- 

0.1818 

Case PlA (Uniform Inlet) 

- 

- 

1.2 

Case P1B 

A = 
an = 

Z 1 = 

0.004 

-0.199 

0.006 

0.3346 

Case PIC 

A = 

an = 

Z 1 = 

0.00001 

0.0445 

0.0015 

0.5536 

Case P2 

A = 

0.055 

2.172 


an = 
B = 
bn = 
C = 

Z 1 = 

z 2 = 

6.0486 
0.05 
-0.6 
3.5924 
0.01445 m 

0.14161 m 




Table 4: Flow Conditions: Radial Flow 


General data: 


Case R1A: 


Case RIB 


Case R2A: 


Case R2B: 


p = 913.026 Kg/m 3 

v = 7.432 x 10' 6 m 2 /s 

r- = 0.0508 m 
m 

r , = 0.1953 m 
out 


h. = 0.02 m, V. = 0.14996 m/s 
in ’in ' 

q - 9.573 x 10~ 4 m 3 /s, g = 0.0 

Free surface profile: h = h^ n ( ) an 


h. = 0.000508 m, V- = 0.1248 m/s 
in ’ m ' 

q = 2.0236 x 10" 5 m 3 /s, g = 0.0 


Free surface profile: 


*> ■ »i. < r f-Ht 

in 

h = C h in [ 2 - (jf-) b “ 

m 


h. = 0.02 m, V- = 0.14996 m/s 
in in ' 

q = 9.573 x 10' 4 m 3 /s, g = 9.81 m/s 2 

Free surface profile: h = h^ n ( r r ~ ) an 


h. = 0.000508 m, V. = 0.1248 m/s 
in ’in ' 

q = 2.0236 x 10' 5 m 3 /s, g = 9.81 m/s 2 

Free surface profile: h = h^ n ( ^ ) an , 


for r < r^ 
] , for r 


for r < Tj 


IV 



Case H.2C: 


Q = 

^out 

Free 


= 0.000508 n, V in = 0.1248 m/s 

2.0236 x 10" 5 m 3 /s, g = 9.81 m/s 2 

= 0.00030265 d (Fr = 1.0) 

surface profile: h = h. ( — — ) an 

m r in - a 

h=Ch in [ 2 - (rT~) b 

m 


h = h out> at r = r 


out 


for r < r^ 

] , for r > r^ 
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Table 5: Summary of Results: Radial Flow 



Problem 

ODtimum Profile 

RMS Error in 
Pressure (Pa) 

Case 

R1A 

A 

_ 

-0.05 

1.238 



an 

= 

-1.0 


Case 

RIB 

A 


0.044 

0.2961 



an 

= 

2.2876 




r l 

= 

0.06236 




C 

— 

9.74 




bn 


0.02 


Case 

R2A 

A 

_ 

0.01 

1.06 



an 

zz 

0.04 


Case 

R2B 

A 

_ 

0.04 

0.5478 



an 


8.8337 




r l 

= 

0.05369 




C 

— 

8.1225 




bn 

= 

0.05 


Case 

R2C 

A 


0.04 

0.7552 



an 

= 

8.822 




r l 

= 

0.05369 




C 

zz 

8.1 




bn 

= 

0.0 
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Fig. 1.1 The coordinate system on a grid cell 
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PI. Plane falling film P2. Plane film flow on horizontal plate 



P3. Plane film flow under zero gravity 
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Rl. Radial film flow under ,p'-. zero gravity 



////// /// /7TVy //////// 7 

R2. Radial 1 ; film flow on horizontal plate 


Fig. 1.2 Flow systems in present investigation 
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FigJ.4 Analytical solution for radial flow over a plate in the presence of gravity 
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Fig.15 Comparison of the solution for different grid sizes 
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Fig.1.6 Film height in a developing falling film flow 
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Fig.l 7 Free surface pressure in falling film flow 



CO CNJ O 

(pd) 3ynss3dd 


58 


0.000 0.005 0.010 0.015 0.020 0.025 0.030 

VERTICAL DISTANCE (m) 



Fig.1.8 Friction coefficient in falling film flow 
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Fig.1.13 Friction coefficient for plane horizontal flow 
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Fig.1.14 Film height for plane flow under zero gravity 
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Fig.1.15 Free surface pressure for plane flow under zero gravify 
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Fig.1.16 Velocity profile for plane flow under zero gravity 
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Rg.1.17 Wall shear stress for plane flow under zero gravity 
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Fig/ 1.18 Friction coefficient for plane flow under zero gravity 
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Fig.|.20 Velocity profile for radial flow (Cases R1A and R2A 

case R1A, z/zm=0.4 case f^A^Z/^m=(X4 

case R1A, z/zm=0. 8 case R2A, z/zm=0.8 



wA/A ‘33NVlSia IVtNdON 

71 


0.00 0.05 0.10 0.15 0.20 



Rg.].21 Free surface pressure for radial flow (Cases R1A and R2A) 
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Fig.1-22 Wall shear stress for radial flow (Cases R1A and R2A) 
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Fig. 124 Film height 
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Fig.l-25 Velocity profile for radial flow (Cases RIB and R2B) 

case RIB, z/zm=0.4 case 

case RIB, z/zm=0. 8 case R2B, z/zm=0.8 







Fig.1.26 Free surface pressure for radial flow (Cases RIB, R2B and R2C) 
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Fig.1.27 Wall shear stress for radial flow (Cases RIB, R2B and R2C) 
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Fig.|.28 Friction coefficient for radial flow (Cases RIB, R2B and R2C) 
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One-dimensional plane flow analysis under 1-g and 0-g 


2.1. Introduction 

The governing equations and boundary conditions are presented for the situation of a 
thin liquid layer emanating from a pressurized vessel and traveling along a horizontal plate 
with a constant initial height and uniform initial velocity as shown in Fig. 2.1. This situation 
is the same as channel flow, but since the liquid height is very thin the effect of viscosity 
must be accounted for. It is desired to find the liquid height at any distance down the length 
of the plate for different Froude numbers and Reynolds numbers specified at the inlet. Since 
the inlet Froude number will usually be greater than unity, it is possible that a hydraulic 
jump will occur at some point in the computational domain. A hydraulic jump is when 
the flow suddenly changes from supercritical (Fr > 1) to subcritical (Fr < 1) flow, which is 
accompanied by a sudden increase in the liquid height. This is analogous to the shock wave 
in gas dynamics when the flow changes from supersonic (M > 1) to subsonic (M < 1) flow 
in a very short distance. The similarity between the hydraulic jump and the shock wave 
in gas dynamics suggests using the familiar approach of modeling the flow as a transient 
phenomena and allowing the solution to march in time to achieve the desired steady state 
results. The effects of microgravity on the flow is also examined. 

2.2. Mathematical Modeling 

The generalized governing equations for an incompressible fluid with constant properties 
are as follows: 

Continuity equation: 


dp ^ d 
dt dx 


( pu ) + + 


d_ 

dz 


( pw ) = 0 
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y-direction momentum equation: 


( dv dv dv dv\ dP 

— — f u — — |-d- — (- in— = — — — 

\dt dx dy dzj dy 


d T iy ^Tyy dT Z y \ 

dx dy dz ) 


+ P9y 


x-direction momentum equation: 


( du du du du\ dP ( dr xx dr xy dr zx \ 

"l at + u ai + v dy ~ w Tz) = “ Tx ~ [~a^ + ~ai + ^Tj + "• 

The general boundary conditions for these equations are those on the flat plate and on 
the free surface of the liquid: 


y— 0: u=0, v=0 


y ~ h: 


y=h: P + 


1 + 


cr dh 

V*72 Si 


ah 


2pu{\ + t> 2 ) du _ a 

li-W 


y = h: Pnr = 0 

where h is the film thickness; p nr is the stress tangent to the film surface; and a is the 


surface tension. 


The following assumptions are introduced to reduce the complexity of the governing 
equations: 

Assumptions: 

• Incompressible fluid 

• w = Wz = 9 * = 0 

• Boundary layer assumptions 

• No surface tension 

• No interfacial shear on the free surface 

The governing equations and boundary conditions then reduce to the following forms: 
Continuity equation: 


81 



y-direction momentum equation: 


du dv 
1 - — = 0 

dx dy 


x-direction momentum equation: 

du du du 1 dP d 2 u 

dt + U dx V dy p dx ^ dy 2 

Boundary conditions: 


P(x, h,t) = 0 
u(xi,y,t) = Ui 
u(x,y, 0) = u 0 



v(x,0, t) = 0 

Integrating the conservation equations across the thin layer results in quasi - one- 
dimensional governing equations: 

Continuity equation: 


du dv 
dx dy 


Leibniz’s rule: 



= 0 


dx 


f f(x,t)dt = / 
J A J A 


B 

A 


df(x,t ) 
dx 


, N dB ,, dA 
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r h (du\ d r h 

jo {ax) dy= rJo 


a r h 


dh 

dx 


The kinematic condition at the free surface for unsteady problems is used as one of the 
free surface boundary conditions: 

dh dh 

m + 


dh 


dh 




a [ h ( dh\ 

BiJo uiy -\ vlh ~ m) + * 1 * “ 1,10 = 0 


u 


= 1 hJo uiy 


d . LN dh n 

+ di = 0 


y-direction momentum equation - 


( 2 . 1 ) 


dP 

* y =~ P9 



P = -pgy + c 


P = pg(h - y) 


( 2 . 2 ) 
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x-direction - momentum equation - 


du du du 1 dP d 2 u 

dt + dx dy p dx + V dy 2 


/:©* * /.'(■£)**/ (-g)* - -;r (£ w ©>* 




/ du 


f h fdu\ d r h 

L [m) dy ~aiL ud y ~ ulh ' 


ah 

at 


«,- M i ah 

Tl (nk) -uW~ 


dt 
d r h 


ah 

dx 


f n ( ap\ a r ft 
Jo [ai) iv= alL 
a r h 

~aiJ 0 rs( h -y)dy 

-£«»*■) 

=/:<g) 


ai (fiA) { u rF y+ i ^r dy -t 


h 

dy) 


From the continuity equation 


l { u ^) dy = -l\ ua £) dy 


a_ 

Ft' 


(nh) - n\ h ■ ~ + J o (u^jiy + («»)|» - (u»)|o + jf (^)dy 
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0 



{uv)\h = u\h • Hh 



= (tt| h ) 2 • 


dh dh 
dx dt 
dh 

w- + u \h- 


) 

dh 

dt 



dy 


d (-h\ I 8h 
Ft (uh)-u U- w 


'T9(u 2 j 


j [ 


dx 


, , , 2 dh dh d ( 1 _ 2 \ / du\ 

iy + M "gi +u ^"m = ~9i\2 sk j 


dt dx\ 2 * 


\dyj 


Io{^ar} dv = rxIl v?iv - ( u| ‘ )1 


dh 

dx 






5i// 


a 


(u/i) + — / u 2 dy H — 51/1 2 

fir } dxlJo y 2 y J 


du\ 

= - U 'Ty) 


(2.3) 


The following assumptions were made to simplify the governing equations: 



w T ith an analogy to the Blasius solution of frictional flow over a flat plate where Moo is replaced 
by u: 


In addition: 



0.332u 2 

sjuxjv 
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dh 
. dt 


d _dh du du 

— (uh ) -- u — — h h — = h — 
dt V dt dt dt 


which resulted in the following governing equations: 


uh = Q = constant 


(2.4) 


du 
dt dx L 


u*h + -an 

2 


0.332u 2 

sjux/u 


(2.5) 


These equations form a system of nonlinear partial differential equations. After the 
equations are nondimensionalized, the dimensionless liquid height is eliminated from the 
momentum equation which leaves one nonlinear partial differential equation with one de- 
pendent variable which is the dimensionless velocity. The assumption involving the term 
dh/dt was made to simplify the governing equations so that the height of the liquid film 
could be eliminated. By using this assumption, it is understood that the unsteady solutions 
with respect to time are not accurate, but the accuracy of the steady state solution is not 
affected. An analysis concerning the term Jq u 2 dy is given in Appendix 2.1. 

Dimensionless variables: 


u 

Uj 


= V 







= Re t 


u = uiV 


h — h\8 x = h\£ 


t = 


h\r 

Ul 


x-direction momentum equation: 
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h(u) t + ( u 2 h + ^ gh 2 ) x = 


0.332u 2 



(M) 


d(uiV) 

d(h}r/u\) 


d 

d{h\i) 


(u,V) J (M) + js (M ) 2 


= -0.332(m]V) 2 


(fiiV)(AiO 


■u. 


c dv 2 




V 2 8 + 


sM 2 \ 


2w 2 


— 0.332U] 



Wr+ (V 2 5 + 



0.332F 2 

V^TW 


Continuity equation: 


0.332V' 3 / 2 


uh — Q 


(■ uiV)(hi6) = Q 

VS = constant =1 or 6 = Substituting and dividing by | gives: 


(V), + irv+ _i L') _ lf -0.332W^ 

1 >, + A + 2Fr\ VV e A v^T l ) 


(V), + v( 


v + 


1 1 


^ = 


0.332V' 5 / 2 


2 Fr\ V* / ( v'Eei? 


Dimensionless governing equation: 


V r + VG i = H 


( 2 . 6 ) 


87 


w 


here 


G = V 


1 ] 
2Fr\ ' V 7 * 


H = 


0.332F 5 / 2 


\/Rei( 

For the case when the gravity is zero, the equation for G is as follows: 


(2.7) 


G = V 

The governing equation has one time derivative and one space derivative, so an initial 
condition and a boundary value are needed. 

Initial condition: 

r = 0 : V = V 0 = 1.0 for 6 < ( < ( n 

Boundary condition: 

( = (i : Fr = J-Vj, V = 1 

2.3. Numerical solution procedure 

Due to the similarity between the hydraulic jump and the shock wave in gas dynamics the 
MacCormack explicit method, which is quite often used for the solution of compressible flow 
problems, will be used in the present numerical analysis. Since it is an explicit method, the 
unknown variables are found in terms of known quantities, as opposed to implicit methods 
which must solve a matrix equation to obtain the solution of the problem. 

The MacCormack method is a two-step scheme that uses first-order finite-difference 
approximations. The scheme first predicts the solution at the next time step with a forward 
time, forward space (FTFS) differencing scheme and then corrects the prediction with a 
forward time, backward space (FTBS) scheme. This results in a approximation that is 
second-order accurate in both time and space (Hankey, 1982). 

When applied to the one- dimensional linear convection equation, 


du du 
dt + ° dx 


= 0 
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the MacCormack method yields (Anderson et al., 1984) 


u" +:1 — u j — c-— (u" + j - u") Forward predictor 


At 


Ax 


,n+ 1 


+ u" +1 


Af 

Ax 


K +1 - 


Backward corrector 


The backward corrector can be seen as an arithmetic average of the old solution, u”. and 
the new solution based on the predicted solution, — c^(u" +1 - u”*j), which uses a 

backward difference. 


The governing equation for the present problem is as follows: 


V T + VG t = H , G = V + 


Forward predictor: FTFS 


1 1 


2 Ft\ V 2 


H 


-( 


0.3321 75 / 2 

y/Reit 


V™ +1 _ 

+ K n 

Ar j 


rC?+i ~ G] 


A£ 


= H ? 


At 


vr = - v? • + at ^> 


T/n+l _ T/»> 
J J 


1 - S? (G ”+’ “ 


+ Ar//; 


(2.8) 


Backward corrector: FTBS 

The finite-difference equation based on the predicted solution using a forward time, 
backward space differencing scheme is as follows 


(l/f+l)' _ yn+1 


At 


V. n+1 

j 


qti - fl _ grn+] 




= 


89 


At 


0, )' - lj"+ 1 - . 1'"+' (c"" _ G"+Ij + ArJ/V+' 


(v>" +i )' = >7 +i [1 - |j(gj +i - G”+>y 


+ At //” + 1 


The corrected solution is the arithmetic average of the past and predicted solutions 


v-~ l = 2 [y " + (F / +1 ) ,] 


V n +i 

J 


\{ v r+ v r' 


- i -^( o r , - o P0 


+ At//; 


? +1 } 


(2.9) 


where 


G n = V n + 1 ( 1 Y 

1 ' 2 Ft]\V?> 

(2.10) 

„„ 0.332(V7) 5 / 2 

(2.11) 

G r = + 2 ; ( L,) 

1 j 

(2.12) 

0.332(Vp)VJ 

' fy 

(2.13) 

Since the forward predictor velocity is in terms of a forward-space approximation, an 


outlet boundary condition on the velocity is needed. For the case of 1-g, it is assumed that 
the Froude number at the outlet is unity, which is a common boundary condition when a 
liquid falls over an edge because the liquid accelerates from a subcritical velocity to the 
critical velocity. The Froude number and the dimensionless velocity are related as follows: 
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- 2 - 9 

r 2 r 2 u 

Fr ' = s -*7 Fr - s 


/V _ (ft! X^) 2 gfej _ y 2 
Fr 2 g{h]6) u 2 6 


V = 


Fr \ 2 / 3 
FrJ 


k = ^v 7 ’ = t± 

" VFrJ \Fr 


2/3 


E , - 2/3 
r rj 


Frj / \Frj . 

For the microgravity case, the slope at the last node is set equal to the slope at the next 
to last node. 


The solution of the governing equation using the MacCormack explicit method proceeds 
as follows: 

• The parameters pertaining to the numerical domain and the inlet and outlet boundary 
conditions are specified. 

• The initial velocity distribution is input to the program. 

• The variables G and H are computed using the velocity profile at the old time step as 
given in eqs. (2.10, 2.11). 

• The velocity distribution at the midpoint time step is calculated in terms of the velocity, 
G, and H at the old time step as shown in eq. (2.8). An outlet boundary condition is 
needed at this step because of the forward-space approximation. 

• The variables G and H are computed again by using the velocity profile at the midpoint 
time step as presented in eqs. (2.12, 2.13). 

• The velocity distribution at the new time step is calculated using eq. (2.9). The inlet 
boundary condition is used in this step because of the backward-space approximation. 
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• The values of the velocity distribution at the new time step are used as the initial velocity 
profile for the next iteration. 

• The process is repeated until steady values are reached. 
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PSEUDO CODE FOR CHANNEL FLOW 


* Input parameters * 

Fr u Rej, V], 6 , £ n , At, A£ 

* Input initial velocity distribution at the old time step * 

VOLD(l) = 1 

* Compute G and H at the old time step * 

10 GOLD(I) = GOLD(VOLD{J), Frj) 

HOLD(I) = HOLD(VOLD{I), Re u ((I)) 

* Compute V at the mid time step * 

VMID(I) = V MID(VOLD(J), At, A(, GOLD{I + 1), H0LD{1 )) 

* Compute G and H at the mid time step * 

GM1D(J) = GM1D(V MID(I), Fr } ) 

HMID(I) = HMID(VMID(I), Re } , £(/)) 

* Compute V at the new time step * 

VNEW{1) = VNEW{V0LD{1), VM1D{1 ), GM1D(J), GMID(I - 1), HMID(I)) 

* Test for convergence * 

* Let V at the new time step be V at the old time step * 

VOLD(I) = VNEW(I) 

GOTO 10 

* Print results * 
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2.4. Analytical Solution for Microgravity 

The gravity term in the conservation of axial momentum equation can be set to zero 
to simulate a microgravity situation. For one-dimensional steady flow, the conservation of 
mass and axial momentum are as follows: 


0 = uh = constant 


where 



cy u* 

TI 


Setting g = 0 results in: 


Cf 


=0.664 

=0.664 



du Cfu 2 

dx 2 h 

The liquid film height can be eliminated from both equations with the following result: 


u 


-2 



Integrating this equation gives: 


£l 

2 Q 


x — 



u 


The integration constant is found with the following boundary condition: 


x — x\ : u = U] , h — h) 
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u 


Ul 



h 

h 


This equation is solved for h: 


A - 4- (] - Cf Xi ^ 

hi 2 hi xi \ 2 hi) 


h = 




hi 



— xi) + hi 


0.332, 


I vh 
xuihi 


(x - xi) -f hi 


This equation has been solved for the following conditions: 

" = >-55 x JO" 3 


= 3.93 


m 
sec 

hi = 0.01 in 


xi — 2.0 in 


which corresponds to an inlet Reynolds number of Rej = 25.36. 

2.5. Results and Discussion 

The results of the computer program for 1-g are shown Figures 2.2 - 2.7. The figures 
show the converged solutions of the problem with inlet Froude numbers ranging from 1.0 to 
6.0 and inlet Reynolds ranging from 12.68 to 76.06. In the 1-g situation, it is observed that 
a hydraulic jump is predicted and the location of the jump changes with the Froude number. 
At very low Froude numbers, the jump is located almost at the liquid inlet. As the Froude 
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number increases, the jump is pushed back away from where the liquid enters the numerical 
domain, which is what has been seen experimentally. 

Figures 2.8 - 2.10 present solutions of the problem under microgravity conditions for 
inlet Reynolds numbers of 25.36, 50.71, and 101.4. Figure 2.8 also presents the analytical 
results of the case under 0-g. For the microgravity simulation, the results are quite different 
than the results of the 1-g simulation because a hydraulic jump is not predicted. The liquid 
film height increases monotonically along the length of the plate. It is shown in the graphs 
that as the Reynolds number increases, the liquid film height becomes smaller. Figure 2.8 
shows the excellent agreement between the analytical and numerical solutions for 0-g. 

Figure 2.11 presents the dimensionless liquid film height for the case when the slope at 
the outlet is set equal to the slope of the node just inside the computational domain. It was 
found that the solution did not converge with this outlet boundary condition. 

Figures 2.12 - 2.15 present the results of the 1-g simulation when the friction at the 
w r all is increased by a factor of three. The graphs show that as the friction increases, two 
phenomena occur: the thickness of the film increases, and the location of the hydraulic jump 
is closer to the entrance than when the friction is less. This is what one would expect in a 
physical situation. 
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Appendix 2.1 


The following is an analysis of the validity' of the assumption concerning the term Jq v?dy. 
The following approximation is given: 


,h 

J u 2 dy ~ u 2 h 

whereas by the definition of the mean velocity: 

a = k [I« uiy f 

It is therefore necessary to determine whether the following equation is accurate: 

fh - , J r rh -j 2 

/o " iy = hlJo uiy . 

Since the flow situation is similar to boundary layer flow, a modified approximate solution 
to the boundary layer velocity profile will be used to examine the approximation. 




u= 

2 h 

2 h 
dy = —dU 

7 r 


£n 2 i,^Vlj o rl \i^u(^iu) 

2hV 2 [*/* , 

^ / sin 2 UdU 

7T JO 


2hVl ri i 

= - -sin2U 

7T L2 4 

hV 2 

= =0.5 hVl 


1 "V 2 
0 
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For this particular velocity profile, the approximation for the term Jq u 2 dy is within 20% 
of the exact value. Therefore, it is recommended that the approximation be used in the 
present analysis. 
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Fig. 2.1 The coordinate system for channel flow 



100 



xi Dimensionless liquid height versus dimensionless distance 

for channel flow with FR1=1.0, RE1=12.68 , = /, 0 
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Dimensionless liquid height versus dimensionless distance 
for channel flow with FR1=3.0, RE1=38.04 , F f/I -- /. a 
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Dimensionless liquid height versus dimensionless distance 
for channel flow for microgravity simulation with RE1=25.36 
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Dimensionless liquid height versus dimensionless distance 
for channel flow for microgravity simulation with RE1=101.4 
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One-dimensional radial flow analysis including solid-body rotation under 1-g and 0-g 

3.1. Introduction 

The governing equations and boundary conditions are presented for the situation of a 
thin liquid layer emanating from a pressurized vessel and traveling radially along a horizontal 
disk with a constant initial height and uniform initial velocity as shown in Fig. 3.1. It is 
desired to find the liquid height at any distance along the radius of the disk for different 
Froude numbers and Reynolds numbers specified at the inlet. Since the inlet Froude number 
will usually be greater than unity, it is possible that a hydraulic jump will occur at some 
point in the computational domain. A hydraulic jump is when the flow suddenly changes 
from supercritical (Fr > 1) to subcritical (Fr < 1) flow, which is accompanied by a sudden 
increase in the height of the liquid. This is analogous to the shock wave in gas dynamics 
when the flow changes from supersonic (M > 1) to subsonic (M < 1) flow in a very short 
distance. The similarity between the hydraulic jump and the shock wave in gas dynamics 

suggests using the familiar approach of modeling the flow as a transient phenomena and 

allowing the solution to march in time to achieve the desired steady state results. The 
effects of microgravity and solid-body rotation on the flow are examined. 

3.2. Mathematical Modeling 

The generalized governing equations for an incompressible fluid with constant properties 
are as follows: 

r-direction momentum equation: 

du du v du v 2 du\ dP / d 2 u 1 du u 1 d 2 u 2 dv d 2 u\ 

dt dr + rdB r ~^ U dzj dr ^ \ dr 2 ^ r dr r 2 ^ r 2 dd 2 r 2 dO ^ dz 2 ) 

0-direction momentum equation: 
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p 


dv dv v dv uv dv\ 

dt 8r r 86 r dz) 


1 dP , ( d 2 v 1 dv 


r 89 


'+p 


dr 2 r dr 



1 8 2 v 2 du d 2 v\ 

^dfi + ^ d9 + Ih 2 ) +P " 


z-direction momentum equation: 


/ dw dw v dw dw\ dP ( d 2 w 1 dw 1 d 2 w d 2 w\ 

p \m + "a? + ; + ” Tz ) = ~ al + H 9^ + ; a7 + ? W + as) + P3 ‘ 

continuity equation: 

du u 1 dv dw 

dr r ^ r 86 + dz 

The general boundary conditions for these equations are those on the disk and on the 
free surface of the liquid: 

z=0: u=0. w=0, v=o>r 
„_u. ,,, _ dh , 8h 

Z-h. W- 

Vnr = 0, p nn = -ph + f 

where h is the film thickness; p„ r , stress tangent to the film surface; p nn , normal stress; a, 
surface tension;. 

The following assumptions are introduced to reduce the complexity of the governing 
equations: 

• Incompressible fluid 

Q 

• = 9 t — ge — 0 

• Boundary layer assumptions 

• No surface tension 

• No interfacial shear on the free surface 
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• Solid body rotation, v = u >r, but neglect the contribution made by the azimuthal mo- 
mentum equation. 

The last assumption was made as a first step to include the effect of rotation on the flow. 
In the future, the azimuthal momentum equation will be included in the analysis. 

The governing equations and boundary conditions then reduce to the following forms: 
continuity equation: 


du u dw 

I 1 = 0 

dr r dz 


z-direction momentum equation: 


dP_ 

dz 


= ~P9 


r-direction momentum equation: 


du du 


du 


1 dP 


d 2 u 


at + u e-r +w Tz~ ru> = +, 7 m 

Boundary conditions: 

P(r, h, t) — 0 
u(ri,z,t) = Ui 
u(r,z, 0) = uo 

(&)L=° 

w(r,0,t) = 0 

Integrating the conservation equations across the thin layer results in quasi 
dimensional governing equations: 

Continuity equation: 


- one- 


117 


Leibniz’s rule: 


du u 
dr ^ r 


dw 

+ 5l = ° 




dz = 0 


d_ 

dx 


fB /• 1 

L !(z ' l)dt = J A 


df(x,t ) 
dx 


dt + f(x,B ) • 


dx 




dA 

dx 



The kinematic condition at the free surface for unsteady problems is used as one of the 
free surface boundary conditions: 

dh dh 

m + “I* ' Yr = ”1* 
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dz — w\h — ta|o 


l (sA) + 


dh uh 

— - w\ h + — + w\h - uj|o = 0 


9 uh dh 

dr v ' r dt 


1 d(ruh ) dh 

r dr + dt 


(3.1) 


z-direction momentum equation: - 


dP 

JI = - pg 



P = -pgz + C 


P = pg(h - z) 


(3.2) 


r-direction momentum equation: 


du 

It 


du 



du , 

+ w— ror 

dz 


1 dP d 2 u 

L y 

p dr dz 2 



1 

P 



dz + u 



dz 
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h ( du 


a 


_d_ fh 
dt) dt 


[ h , dh 

Jo udz - ulh 'dt 


d. dh 


t h ( du\ r h d(uw) r h ( dw\ , 

Jo \ W Tzr = Jo ~dz~ iz - A 


fh f dw\ 

~(uw)\h - {uw)\o - J \u— jdz 

=(“”)''■ - /„ (“S)^ 

f h /dP\ d f h dh 

L l pdz ~ p ' *-Tr 

d f h 

= Jo P9(h - z)d 2 


r h /d 2 u'' 

)dz= [ 

Jo V dz 2 j 

' Jo 

-Jo 


'du'' 


“(s) 
--(£) 


/du' 
h \ dz ) 


d_ 

dt 


(uh)-u\i, '~gl + j 0 (u£)iz+(uw)\ h - (u-£}d,-rhJ = 
From the continuity equation: 


du\ 
dz) o 


- fo{ u lh) iz = Io( u %) dz + IXt )* 2 

Using the kinematic condition at the free surface: 


, ,, , (dh dh\ 

{uw)\ h = u \h[Qi+u\ h —) 

dh . 2 dh 

=u\ h ■ Tt + wo - Tr 


120 




121 



Assumption: 


rh 

I u 2 dz ~ u 2 h 

Jo 


Jt [uh) + - "(S) 


13 

r 3r 


( ru 2 h ) = - 


5 

r-—{u 2 h) + u 2 h 
or 


a_ 

at' 


(uh) + ~(i 2 t) + | (^A 2 ) + ^ - <-ftu , 2 = 


5 , n 2 ft. 

: s (uA > + — 


u 2 /i 


-(=) 


_a 

at' 


(M) + ^(“ 2/l + ^ 2 ) 


ti ru; 2 — 


n 2 


-'(s)L 


Assumption: 


u no = n 


au\| 0.332n 2 


a 

at' 


d (_ 2 1 , 2 \ / 2 u*\ 0.332 n 2 

^ + a^+2^ j = H w 


u 2 \ 0.332 u 2 


Assumption: 


dh 

at 


2r 0 


a , a/i an an 

—(uh) = u— + h — = h — 

ar J at at at 
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r-direction momentum equation: 


h Tt + rr { a 2 h + \ 9 h 2 ) = h { TU , 2 -' v )- 


Continuity equation: 


1 d . _ , . 
-—(ruh) = 0 
r or 


0.332H 2 

yjurfv 


(3.3) 


ruh = constant = Q 


(3.4) 


Dimensionless variables: 


u 


Ul 


V 


hi 


r ui 

hi = t % = T 


hi o 
tv— = it 

Ul 


Uyhi 


Re i 


v 


u 


gh Fr * 


u = u i V h = h\8 r — h] £ t = tv — 

v,\ h\ 


Radial momentum equation: 


‘I l + TX* 2h + \ 9h2 ) =h (™ 2 -T) 


u 2 \ 0.332u 2 


yj ur ju 


, t CS 9{uiV) , d 

\h\b) 7777 ' ~r~ + 


= (hj8) 


d{h\T ju\) d(hi£) 


(uiF) 2 (M)+^(M) 2 


0.332(uiC) 2 / ^^^ 
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2 \c9V , _ 2 \ d 

( “ l){ a/ + lu > >8( 


vn + l ^ s 1 

2 ui 


- (*?)« 


1/2 


£D 2 - 


(ii] )0.332V /Z / jRe\V^ 


s ir + §*\ v2s 

OT oi L 
Continuity equation: 


2Fr\ 


= 6 


l r 2 1 


£D 2 - y j - 0.332 vW/jRat 


ruh — O 


(hiOfaVHhjf) = Q 


Q 




Divide both sides by 6 and eliminate it from the radial momentum equation with the 
continuity equation: 


ar i_a 
dr + s at 


v 2 s + 


2 Frf J 


= w 2 


V 2 0.332V 3 / 2 

1 


6 = 


V( 


dV Vi d 


U#(b\, 1 (t i\ 2 1_,o2 V 2 (Vi\0.332V 3 / 2 

r W + 2friWr p ‘ £ U « 


t V 6 / 


av » r . a 

ar + v ^a£ 


e 2^2 J 


V 2 0.332V 3 / 2 ji 

1 6 


Dimensionless governing equation: 
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V, + V(G £ : H 


(3.5) 


where 


G = 


— + 

( 2 Fr\ \ V(J 


(3.6) 


H = £fi 2 - 


0.332F ' 


f 6 V *ei 

For the case when gravity is zero, the equation for G is as follows: 


(3.7) 



The governing equation has one time derivative and one space derivative, so an initial 
condition and a boundary value are needed. 

Initial condition: 


r = 0 : V — Vo = 1.0 for < £ < £ n 

Boundary condition: 

£ = £i : Ft = Fr\ , V = 1 

3.3. Numerical solution procedure 

The MacCormack explicit method is quite often used for the solution of compressible 
flow problems. Since it is an explicit method, the unknown variables are found in terms of 
known quantities, as opposed to implicit methods which must solve a matrix equation to 
obtain the solution of the problem. 

The MacCormack method is a two-step scheme that uses first-order finite-difference 
approximations. The scheme first predicts the solution at the next time step with a forward 
time, forward space (FTFS) differencing scheme and then corrects the prediction with a 
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forward time, backward space (FTBS) scheme. This results in a approximation that is 
second-order accurate in both time and space (Hankey, 1982). 

When applied to the one-dimensional linear convection equation, 


du du 

at +c ai = 0 


the MacCormack method yields (Anderson et al., 1984) 


,n+l 


At 


uj ' ' = u" — c— -(u” +1 — u") Forward predictor 

j" 


.n-fl __ 


1 


At 

Ax 


K +1 - U 1- 1) 


Backward corrector 


The corrector can be seen as an arithmetic average of the old solution, u", and the new r 


solution based on the predicted solution, u” +1 — — it"**), which uses a backward 

difference. 


The governing equation for the present problem is as follows: 


V T + V(G( = H 


■where 


G = 


^ + Jl_MV 

£ 2Fr\ \V£) 


H = W 2 


V 2 

T 


0.332V 5 / 2 /~£~ 

6 


Forward predictor: FTFS 


i/n-t -1 _ T/n 

7 J 


Ar 


+ 'TO 


G n 

A£ 
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( 3 . 8 ) 


A r 


V , ," +I = '? - (G?J + , - CJ) + At//” 


i/n-fl _ T/n 
J > 


Ar 


1 - 6-^(G>+l - G?)J + Arif" 

Backward corrector: FTBS 

The finite-difference equation based on the predicted solution using a forward time, 
backward space differencing scheme is as follows 


v rugn + ^ 


G n + 1_ G „ + 1 


Jbl 


A£ 




( vj * +i )' = v; +i - vr^^r 1 - G i~i) + An/- 1 


(^ n+1 )' = l T +1 [l - 6 ^ (cf 1 - G ”-i )] + At i /" +1 

The corrected solution is the arithmetic average of the past and predicted solutions 


V ? +1 = ±[V? + {vr+'fl 


V ?* 1 = + V ?* 1 


1 - 6 fi( G T'- G ”- + . 1 )] +Ar ' ff i' +1 } 


where 


( 3 . 9 ) 


r* V i , fa / 1 V 
' 0 2 Fr|W/f,J 


( 3 . 10 ) 


H] = (jll 2 - 


(V /*) 2 0 . 332 (lj") 5 / ! rj~ 


0 


6 V Bei 


( 3 . 11 ) 
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(3.12) 


V '" +1 

G n + 1 _ j 


+ 


fc /_L-Y 

2 Fr\ \ 




= 


(> 2Fr2 


(V" +1 ) 2 0.332(V / / +1 ) 5 / 2 rjf 

tj 6 V ^ e i 


(3.13) 


Since the forward predictor velocity is in terms of a forward-space approximation, an 
outlet boundary condition on the velocity is needed. For the case of 1-g, it is assumed that 
the Froude number at the outlet is unity, which is a common boundary condition when a 
liquid falls over an edge. The Froude number and the dimensionless velocity are related as 
follows: 


Fr 


Ft 2 = - 


gh\ 


.-.2 


gh 


Fr 2 (nj V) 2 ghi V 2 


Fr 2 


g(hi6) 


u; 




Ay 3 

6 


Mi-sr 




1/3 


For the microgravity case, the slope at the last node is set equal to the slope at the next 
to last node since the outlet boundary condition can not be specified by the Froude number. 

The solution of the governing equation using the MacCormack explicit method proceeds 
as follows: 

• The parameters pertaining to the numerical domain and the inlet and outlet boundary 
conditions are specified. 

• The initial velocity distribution is input to the program. 
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• The variables G and H are computed using the velocity profile at the old time step as 
given in eqs. (3.10, 3.11). 

• The velocity distribution at the midpoint time step is calculated in terms of the velocity, 
G, and H at the old time step as shown in eq. (3.8). An outlet boundary condition is 
needed at this step because of the forward-space approximation. 

• The variables G and H are computed again by using the velocity profile at the midpoint 
time step as presented in eqs. (3.12, 3.13). 

• The velocity distribution at the new time step is calculated using eq. (3.9). The inlet 
boundary condition is used in this step because of the backward-space approximation. 

• The values of the velocity distribution at the new time step are used as the initial velocity 
profile for the next iteration. 

• The process is repeated until steady values are reached. 
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PSEUDO-CODE FOR RADIAL FLOW 


* Input parameters * 

Fr 1 , Re 1 , £ 1 , ( n , A r, A£, 

* Input initial velocity distribution at the old time step * 

VOLD(I) = 1 

* Compute G and H at the old time step * 

10 GOLD(I) = GOLD{V OLD(l), Frj) 

HOLD(I) = HOLD(VOLD(I), Re u £(/)) 

* Compute V at the mid time step + 

VMID(I) = VMID(VOLD(I), At, A{, GOLD(I + l), H0LD{1 )) 

* Compute G and H at the mid time step * 

GMID(I) = GMID(V MID(I), Ft j) 

HMID(I) = H MI D(V MI D(I), Re u £( I )) 

* Compute V at the new time step * 

VNEW(I) = VN EW(VOLD(I), VMID{1), GMID(I), GM1D{1 - 1), HMID{I )) 

* Test for convergence * 

* Let V at the new time step be V at the old time step * 

VOLD(I) = VNEW(I) 

GOTO 10 

* Print results * 
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3 • V Results 


The results of the computer program that predicts the height of liquid in radial flow 
with a free surface on a disk under 1-g are shown in Figures 3.2 - 3.6. These figures show 
the dimensionless liquid film height under 1-g with no rotation as the initial Froude number 
changes from Frj = 2.0 - 6.0 and the initial Reynolds number changes from Rej = 25.36 - 
76.06. The outlet boundary condition for these cases is Fr n = 1.0, which corresponds to a 
liquid falling over an edge. For small initial Froude numbers (Frj — > 1), a radial hydraulic 
jump is predicted very close to the inner radius of the disk. A hydraulic jump is when the 
liquid film height suddenly increases when the flow changes from supercritical (Fr > 1) to 
subcritical flow (Fr < 1). As the initial Froude number increases, the radius of the hydraulic 
jump increases, which agrees with experimental data. For all of the cases in this set of figures, 
the liquid film height first decreases and then increases in the supercritical region. As will 
be shown in other cases, this is not always true if the inlet Reynolds number is sufficiently 
small. 

Figures 3.7 - 3.11 give the results of when the gravity term is set to zero, which cor- 
responds to a microgravity simulation. The inlet condition was specified with the inlet 
Reynolds number which ranged from Rej = 25.36 - 76.06. The outlet boundary condition 
was specified by forcing the slope of the curve at the outlet be equal to the slope just before 
the outlet. In all of the cases studied, a hydraulic jump did not appear. The trend in all 
of the cases in these figures was an initial decrease in the liquid height followed by a mono- 
tonic increase of the liquid film height. It can be seen that as the initial Reynolds number 
increases, the liquid film height decreases at all points along the radius. 

Figures 3.12 - 3.14 demonstrate the effect of solid-body rotation on the liquid film height. 
Figure 3.12 shows the flow field without rotation, Fig. 3.13 presents the liquid height with 
the dimensionless angular velocity set to ft — 10 -3 , and the Fig. 3.14 shows the film height 
with ft = 1CT 2 . With ft = 10 -3 , the shape of the hydraulic jump has been ’stretched’ along 
the radial axis, while with ft = 10~ 2 the hydraulic jump has been completely ’washed out.’ 
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Figure 3.15 shows the results of the program when there is gravity, no rotation, and the 
outlet boundary condition is specified such that the slope at the outlet is equal to the slope 
just before the outlet. It can be seen that the solution does not converge even when the time 
step is decreased to 0.1 of the time step with the outlet boundary condition specified. 

Figure 3.16 compares the results of the channel flow with that of the radial flow with 
the same inlet and outlet conditions. The height of the liquid film for channel flow is 
approximately twice that of radial flow. 

Figures 3.17 - 3.19 compare the results of the computer program under 1-g when the 
friction at the wall is increased by a factor of three. The graphs show that the effect of 
increasing the friction is the same as in channel flow: the height of the liquid film is increased, 
and the location of the hydraulic jump is closer to the inlet of the computational domain. 
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